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Abstract 

1//" noises are ubiquitous and affect many measurements. Tliese noises are both 
a nuisance and a peculiarity of several physical systems; in dielectrics, glasses and 
networked liquids it is very common to study this noise to gather useful informa- 
tion. Sometimes it happens that the noise has a power-law shape only in a certain 
frequency range, and contains other important features, that are however difficult 
to study because simple fits often fail. Here I propose a model-based fit procedure 
that performs well on spectra obtained in a molecular dynamics simulation. 



1 Introduction 



1//" noises are very common and affect many measurements; the literature 
on tliis subject keeps growing and the apparent ubiquity of these noises has 
always drawn a great deal of attention. In the experimental practice, they are 
both a nuisance and a pecuharity of several physical systems; in dielectrics, 
glasses and networked liquids it is very common to study these noises to gather 
useful information [1,2,3,4]. Sometimes it happens that the noise has a power- 
law shape only in a certain frequency range which spans several decades, and 
at the same time contains other important features, that are however difficult 
to study because simple fits often fail. The main reason of this failure is that 
the prominent low-frequency peak biases the fit so much that the minute and 
mostly high-frequency features are neglected. Here I propose a model-based 
fit procedure that bypasses this problem and that performs well on spectra 
obtained in a molecular dynamics simulation of water. 
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In the rest of this introduction I review the classic superposition argument 
that relates power-law spectra to the single exponential relaxation processes; 
in section 2 I analyze the properties of the autocorrelation function of 1//°" 
spectra; in section 3 I consider the spectral behavior associated to some well- 
defined distributions of relaxation rates; finally in the last section I show the 
results of a model-based fit in the case of a molecular dynamics simulation of 
liquid water, and I summarize my conclusions. 

A mathematical mechanism for producing 1//" noise was proposed long ago 
by Bernamont [5], who observed that the superposition of many Lorentzian 
spectra with a certain distribution of different rates could produce a spec- 
tral density with a 1// region. The Bernamont superposition argument can 
be made rigorous with a slight modification of the standard proof of Camp- 
bell's theorem [6], and it goes as follows. Take a signal x{t) originated by the 
linear superposition of many random pulses, i.e., pulses that are random in 
time and can be described by a memoryless process with a Poisson distribu- 
tion, have random amplitude A drawn from a distribution with finite variance 
and probability density g^iA), and such that their pulse response function 
h{t,X) = exp(— At) (if t > 0, otherwise h(t,X) = 0) is drawn from a distri- 
bution with probability density gx{X). The pulses are received and detected 
with a rate n{A, A) which in general depends both on the amplitude A and on 
the decay rate A. The pulse arrival process is Poissonian and thus one detects 
on average [n{A, X)dAdX\ dt pulses in the time interval (t', t' + dt) (and in the 
amplitude- A range dAdX)] for the same reason the variance of the number of 
detected pulses is also equal to [n{A, X)dAdX] dt. This means that the mean 
square fiuctuation of the output signal at time t is given by the integral 

^max -^max t 

{{Axf)^ J gx{X)dX J gA{A)dA J dt'n{A,X)[Ah{t-t',X)f (1) 



If we assume that the rate of occurrence n does not depend on A and A, and 
rearrange the time integration, then the integral (1) simplifies to 



{{Axf)^n{A^) J g,{X)dxJdt[h{t,X)f (2) 



Now let H{uj,X) be the Fourier transform of h(t,X), then from the causality 
constraint on h{t, A) and Parseval's theorem we find that the mean square 
fiuctuation (2) can be trasformed into 
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j du j g^{X)dX\H{u,X)f (3) 



The right-hand expression in equation (3) shows that the spectral density is 

A- 



5M = ^ / 9x{X)dX\H{u;,X)\' (4) 



and since \H{laj,X)\ = {u"^ + A^) ^ we obtain eventually 



^^^^ = J ZpTx^"^^ 
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If we assume that the decay rates A are uniformly distributed between X^in 
and Xmax (i-e., gx{X) = {X^ax - Amm)~^ ) the spectral density becomes 

(-1/ \ n(^A ) 1 / X^nx , Xfjiin\ /„\ 
o[u!) — — — -— arctan arctan (6) 

27r(Amaa; - Xmin) \ <^ ^ J 



so that S{u}) is approximately constant if < a; -C A^m -C Xmaxi and it is 
approximately equal to 

n{A') 1 



'^'^i.Xmax XjYiiji) 



if ^min Amax "C u! , whilc it is approximately equal to 
njA') 1 

in the region in between the extreme rates (A^i„ <^ a; <^ A,, 



(8) 



The spectral density (6) has an intermediate region with a 1// behavior, 
however most observed spectra are not quite 1// but rather 1//" with a 

ranging from about 0.5 to nearly 2: how can we obtain such spectra using a 
superposition as above, i.e., sampling a distribution of relaxation processes? 
We could take, e.g., a nonuniform distribution of relaxation processes like 
gx oc A"'^, then in the region A^j„ <^ cu <^ X^ax we would find 
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We shall return to these distributions in section 3. 



2 The rate distribution from the correlation function 



We see that from a given rate distribution wc obtain a certain spectral density: 
can we do the reverse and obtain the rate distribution from a given spectral 
density? This is not obvious because the spectral density is only a second- 
order statistics, and does not contain phase information (nor is it possible to 
preserve it for a noise process). However the answer is yes, the rate distribution 
can be recovered from the spectral density. This can easily be seen from the 
formal Taylor expansion of the denominator in the integral (5): 




This expansion is only formal inasmuch as it does not converge everywhere, 
however it shows unequivocally that the shape of S(uS) depends only on the 
even moments about the origin of the probability density g\. A probability 
density function is uniquely determined by the knowledge of all the moments 
(A"^) (see, e.g., [7]), and the even moments alone are not enough, but we could 
still do without the odd moments if the probability density function were an 
even function. This is not so, because the decay rates A must be non-negative, 
and thus the associated probability density function does not have any defi- 
nite parity. However a probability density function which is non-zero only for 
positive values of the decay rates can be written in a unique way as the sum of 
an even and an odd function 5'a(A) = gx"^'^\X) + g\^^^\X), where g^x^^\X) = 
'^"^A) = g,{X)/2 if A > and ^^(A) = ^"^A) = -g,{-X)/2 if 
A < 0, therefore the odd moments can be computed from the even moments 
of the distribution, and the even moments alone uniquely identify the rate 
distribution. 

The previous result is only formal and does not yield a practical inversion 
formula; the actual inversion can be performed in the time domain when we 
recall that the spectral density S{u) is related to the correlation function R(t) 
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by the Wiener-Kintchine theorem 
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then we see from equation (13) that the correlation function is also the Laplace 
transform of g\[X)j(2X)^ and the rate distribution function is uniquely deter- 
mined by the spectral density and can be retrieved by means of a numerical 
inverse Laplace transform. In practice, rather than a numerical evaluation of 
the inverse Laplace transform, one is forced to fit a discrete set of decaying 
exponentials, and moreover from the correspondence between the Bromwich 
inversion integral and the inverse Fourier transform, and from the sampling 
theorem, we see that we must sample the time correlation function, and there- 
fore the noise signal, at a frequency at least twice as high as \m,ax to retrieve 
g\. Notice that because of the A in the denominator of the integrand in (13), 
the slow relaxations are more heavily weighted in the integral, and the high- 
frequency parts of the decay rate distribution are much harder to recover 
than the low-frequency parts; this makes even harder an inversion task which 
is already known to be very difficult [8]. 

The mixtures of decaying exponentials that characterize many experimental 
measurements differ significantly only at very short times, while for longer 
times all the exponentials are equally buried in noise. Disentangling the mix- 
ture and finding the relative weights of the different components is possible 
only if sampling times are very closely spaced at the beginning (and one com- 
mon strategy is to use logarithmically spaced sampling times (see, e.g. [9])) 
and only if one includes some form of prior or assumed knowledge of the dis- 
tribution of decay rates. There are very few well-established procedures to do 
this, and the best known are the programs CONTIN and UPEN. CONTIN 
[10] uses the following strategies: a) it takes into account absolute prior knowl- 
edge, i.e. whichever exact information that may be available at the beginning, 
like the non-negativity of decay rates; b) it assumes some statistical prior 
knowledge as well, which is essentially the knowledge of the statistics of the 
measurement noise; c) it uses a principle of parsimony, which is similar to the 
principle of maximum entropy, though not as well defined. UPEN (Uniform 
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PENalty) [11] assumes instead a priori that the distribution of decay rates is a 
continuous function and penahzes distributions which are either discontinuous 
or have wildly varying curvature. 

In addition to constraints on the shape of the distribution function it is com- 
mon to use some well-defined standard functions that appear to fit very well 
many sets of experimental data; the Kohlrausch- Williams- Watts function de- 
scribes stretched exponentials and works well for relaxations in the time do- 
main and similarly the Havrilijak-Negami (HN) function provides good fits 
to spectral data. These empirical functions are well-known, and in particular 
from the HN spectral shape it is possible to compute analytically the distribu- 
tion of relaxation rates [12]. However, even though these functions often give 
satisfactory fits, it would be much better to connect data from experiments or 
numerical simulations to some well-defined, simple distribution of relaxation 
rates, just like the spectral density in equation (6) can be directly related to 
a fiat distribution of relaxation rates: in the following section I give a list of 
such spectral shapes. 



3 A gallery of spectral densities 



The spectral density in equation (6) produces an intermediate region with 
a 1/f behavior, and includes both a minimum and a maximum relaxation 
rate: at a frequency lower than the minimum relaxation rate the spectral 
density whitens and becomes nearly flat, while at a frequency higher than the 
maximum relaxation rate the spectral density bends downward and assumes a 
1//^ behavior, and for fitting purposes we define the standard spectral density 

'S'flat('^; Amin, >^max) = " ( arctau ^^^^ - arctau ^i:^ ] (14) 



However either the minimum or the maximum relaxation rate (or both) may 

be out of the experimental or numerical simulation range: in these cases the 
bends at low- and high-frequency become invisible, and a fit with the spectral 
density (14) is unstable (at least one of the range parameters is invisible and 
the chi-square hypersurface flattens out in that direction, adversely influencing 
the fit) . This can be corrected using the modified spectral density 



arctan 
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(15) 



when the maximum observable frequency is smaller than the maximum relax- 
ation rate (and the minimum relaxation rate is in the observable range). We 
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Fig. 1. Plot of the spectral density (6) (solid line); the dotted, dashed-dotted, 
and dashed lines represent respectively 1//, 1//^'^, and 1//^ spectra. Both spec- 
tral values and frequencies are given in arbitrary units; here Xmin = l{a-u.) and 

Amao: = 1000(0.^.). 
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Fig. 2. Plot of the spectral density (15) (solid line); the dotted, dashed-dotted, and 
dashed lines represent respectively 1//, 1//^'^, and 1//^ spectra. Both spectral 
values and frequencies are given in arbitrary units; here Xmin = l(a-it-)- 

should use instead the spectral density 

'S'flat,B(^; >^max) = - arctan ( I (16) 

00 \ 00 J 



when the minimum observable frequency is higher than the minimum relax- 
ation rate (and the maximum observable rate is in the observable range), and 
finally the spectral density 

Sioyeriii^) = — (17) 



when both the minimum and the maximum relaxation rates are out of range; 
the spectral densities (6), (15), and (16) are shown in figures 1 to 4. Using 
(15), (16) or (17) improves the fit stability but means that the final description 
of the relaxation rate distribution is incomplete. 
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Fig. 3. Plot of the spectral density (16) (solid line); the dotted, dashed-dotted, and 
dashed lines represent respectively 1//, 1//^'^, and 1//^ spectra. Both spectral 
values and frequencies are given in arbitrary units; here Xmax — 

1000(a.u.). 

We have already given a simple argument that shows that a nonuniform dis- 
tribution of relaxation processes like gx oc A'"''^ between the maximum and 
minimum relaxation rates Xmim A„,„;^ . produces a spectral density with an 
intermediate 1//^+^ region: an exact integration yields the spectral density 



(1 — p)uj'^ 



2 ' a;2 




(18) 



where F{a,b;c;z) = Sfelo ■^^j^y^fr hypergeometric function and (3 e 

(— 1, 1). Just as in the 1// case either the minimum or the maximum relaxation 
rate (or both) may be out of the experimental or numerical simulation range 
and a fit with the spectral density (18) becomes unstable, and this can be 
corrected with the modified spectral densities 



S'p1,a(<:^; Xmin, = L{uj, (5) - 
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when the maximum observable frequency is smaller than the maximum relax- 
ation rate (and the minimum relaxation rate is in the observable range) and 
where the function 

v^-P /i — /9 ^ — B — \ 



is shown in figure 4 and is well approximated by the rational function 

TT 1 



(21) 
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Fig. 4. Graph of the function 
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L{uj,f]) (this product depends on /? alone); 



the dots are obtained from numerical estimates of the r.h.s. of equation (20). 
with 



(1) C2 « -1.2337 

(2) C4Ri 0.253669 

(3) C6 ^ -0.0208621 

(4) cg ^ 0.000917057 

(5) cio ^ -0.0000235759 

The spectral density 

5„,,B(^;A„/?)=j^F(i^.l:i^;Z^) (22) 



works when the minimum observable frequency is higher than the minimum 
relaxation rate (and the maximum observable rate is in the observable range), 
and finally the spectral density 

(23) 



when both the minimum and the maximum relaxation rates are out of range 
(here I extend the notation of definition (17) ); the spectral densities (18), 
(19), and (22) are shown in figures 5 to 7. 

In addition to these distributions, it is possible to consider other shapes like 
g\{X) (X a + bX so that the resulting spectral density from equation (5) is the 
sum of a spectral density like the one in equation (6) plus a term proportional 
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Fig. 5. Plot of the spectral density (18) (solid line); the dotted, dashed-dotted, 
and dashed lines represent respectively 1//, 1//^'^, and 1/ p spectra. Both spec- 
tral values and frequencies arc given in arbitrary units; here \min = l(a.u.), 
^max = 1000(a.u.), and f3 = 0.5. 




c.cccc: ■ 



1 M 100 1000 10000 

frequCTcy (a.u.) 

Fig. 6. Plot of the spectral density (19) (solid line); the dotted, dashed-dotted, 
and dashed lines represent respectively 1//, 1//^'^, and 1//^ spectra. Both spectral 
values and frequencies are given in arbitrary units; here Xmin = l{a-u.), and P = 0.5. 
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Fig. 7. Plot of the spectral density (22) (solid line); the dotted, dashed-dotted, and 
dashed lines represent respectively 1//, and 1//^ spectra. Both spectral 

values and frequencies are given in arbitrary units; here X-max — 

1000(a.tt.), and 

(3 = 0.5. 

to 
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Fig. 8. Plot of the spectral density (25) (solid line); the dotted, line represents 
a simple resonance. Both spectral values and frequencies are given in arbitrary 
units; here uq = 100{a.u.), Xmin = 10{a.u.) and X^ax = 50(a.u.), while the simple 
resonance has loq = 100(0.^.) and A = 10{a.u.). 

but I shall not consider them here, since these shapes seem to be far less 
common than the cases discussed above. 



The integral (5) is a sum of functions that decrease for positive, increasing u 
and therefore cannot be an increasing function and therefore no distribution 
of relaxation rates can possibly describe bumps and other small local features 
such as those that arc observed in the spectral densities of glassy systems. 
These features can be described by resonances or by groups of close resonances; 
the simplest choices are a) fixed resonance frequency and flat distribution of 
relaxation rates; b) fixed relaxation rate and flat distribution of resonance 
frequencies. In the case of a flat distribution of relaxation rates between the 
maximum and minimum rates Xmin, X^ax we find 



Sfr{uj; Xmini Xmaxi ^o) = j 



dX 



UJ — UJq 



X^ + {uj- uo)^ 
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(25) 



and similarly in the case of a flat distribution of resonance frequencies between 
the maximum and minimum frequencies uj^im ^max we flnd the spectral den- 
sities (25) and (26) are shown in flgures 8 and 9. 



dUQ 

x^ + iu- ujQy 



1 

A 



arctan 



U — LJ„ 



X 



— arctan 



U — UJr, 



X 



(26) 
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Fig. 9. Plot of the spectral density (26) (solid line); the dotted, line represents a 
simple resonance. Both spectral values and frequencies are given in arbitrary units; 
here A = 10{a.u.), uJmiu = 50(a.n.) andumax = 150(a.u.), while the simple resonance 
has ojQ = lW{a.u.) and A = lQ{a.u.). 

4 Model-based fit of a simulated spectral density 

When fitting spectra it is important to include the variance of spectral data: if 
Sk is the spectral estimate at the A;-th frequency, and if the time-domain data 
are affected by Gaussian white noise, then the spectral estimate of the white 
noise background has standard deviation [13]; this estimate of the standard 
deviation is usually assumed for simplicity even when there are deterministic 
components or the noise is not white. Moreover if the final spectral density 
is the average of M uncorrelated spectra, then the estimate of the standard 
deviation at the k-th frequency is Sk/\fM- I wish to stress that this treatment 
of the spectral variance is only approximate in the case of colored noises, 
but it is assumed nonetheless, because of the complexity of a calculation that 
includes the correlation between different samples in the time domain (see. 



I have tested the simple model-derived spectral densities described in section 
3 on data kindly provided by C. Chakravarty and A. Mudi [15]: the original 
spectral data are shown in figure 10 and correspond to the 230 K curve in fig- 
ure la of reference [16] (see also [17,18,19] for full simulation details). At very 
low frequency the spectrum is rather steep: a simple fit of the low-frequency 
data shows a 1//^ behavior, and thus we can surmise that this is just the 
high-frequency tail of a simple relaxation with a very low relaxation constant 
(this accounts for 2 fit parameters: amplitude and relaxation rate). At higher 
frequency the slope is smaller and Mudi and Chakravarty estimate a spectral 
index slightly higher than 1 [16]: since there is no hint of a downward bend, 
I exclude the full spectral shape (18) and also the reduced form (22), and I 
choose (19) instead, i.e. I include the possibility of a low-frequency fiattening, 
made invisible by the high-frequency tail of the simple relaxation (this adds 
three more parameters to the fit: an amplitude, a minimum relaxation rate, 
and a spectral index /3). The high-frequency bump resembles rather closely 



e.g. [14]). 
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Fig. 10. Spectral data from [16] (230 K data in figure la): mean square fluctuation 
of potential energy vs. frequency. The overall shape is close to a 1//" spectrum, 
but notice the low-frequency steepening of the spectrum and the pair of bumps: the 
low-frequency steepening can be associated to a strong single relaxation, while the 
bumps should correspond to two resonance distributions like those in equation (26). 

the shape in figure 9, and thus it is reasonable to assume that both the low- 
frequency and the high-frequency bumps correspond to fiat superpositions of 
resonances like in equation (26) (each bump accounts for 4 more parameters: 
an amplitude, a relaxation rate, a minimum and a maximum resonance fre- 
quency, but the relaxation rate is taken to be the same in both bumps). The 
resulting 12 parameter model is: 



A34) + alSf^{uj A34) (27) 

Notice that the assumptions on the relaxation rate distributions help keep the 
number of fit parameters rather low. If we tried to fit with a superposition 
of N simple relaxations we would have 2N parameters (one amplitude plus 
one relaxation rate for each relaxation component): with 12 parameters we 
could fit only 6 simple relaxation components, therefore the assumed shapes 
(that correspond to given distributions of relaxation rates and resonance fre- 
quencies) allow for a much more economical fit procedure. In this case the 
spectral data are averages of M = 448 spectra; table 1 fists the fit parameters 
to the data [15] obtained with a standard Levenberg-Marquardt chi-square 
minimization procedure, and figure 11 compares the fit with the data (the a 
amplitude values in the table are in the spectral amplitude units of fig. 10, 
the A's and the a;'s are in cm"-*^, and /? is dimensionless) . 

The model (27) is a function of both relaxation rate and resonance frequency 
and should thus be described by a two-parameter distribution g{X,ujQ) rather 
than g\{X), however if we concentrate on the projection on the A axis, then we 
can consider only the first two terms: the (reduced) A distribution is shown in 
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Fig. 11. Fit to the spectral data from [16] shown in figure 10 (thick black curve). 
The data are shown in light gray in the background; the dotted curves a, b, c, and 
d represent respectively the first, second, third and fourth term of the model (27). 



Table 1 

Fit parameters for the model in equation (27) to the data from [16] 



ai 


11.435 ± 0.880 


Ai 


0.144 ±0.012 


a2 


1.351 ± 0.0059 


Amm,2 


5.226 ± 0.099 




0.327 ±0.013 


as 


0.102 ± 0.002 


^min,3 


32.1 ± 1.0 




64.4 ± 0.5 


04 


0.02534 ± 0.00002 




421.0 ±0.2 




947.2 ±0.1 


A34 


22.77 ±0.11 



figure 12, and is the sum of a delta-function plus an (unbounded) continuous 
distribution. Notice that such a distribution is quite challenging for other fit 
methods, like those implemented by CONTIN and UPEN. 
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Fig. 12. Projection of the two-parameter distribution g{\,u!Q) that describes the 
model (27) on the A axis. The single relaxation corresponds to a delta-function 
(arrow on the left). 

5 Conclusion 

In this paper I have described a model-based fit of power-law-fike spectral 
densities. Like other similar methods, it embodies a priori information on the 

shape of the distribution, but unlike the other methods, the shape is physically 
motivated, and the fits can be efficiently performed with a reduced number of 
parameters. 
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